function V = multislicetovol(rawfile, volfile, xpix, ypix, xtile, ytile, nslices, xres, yres, zres, flip, dtype, odtype)
% multislicetovol  - convert a multislice MxNx1 image file to a XxYxZ 3-D volume
%
% FORMAT:       multislicetovol(rawfile, volfile, ...
%                               xpix, ypix, xtile, ytile, nslices, ...
%                               xres, yres, zres, dtype, odtype)
%
% Input fields:
%
%       rawfile     filename of multislice raw image data
%       volfile     filename of volume to write
%       xpix        X resolution of multislice matrix
%       ypix        Y resolution of multislice matrix
%       xtile       number of slice tiles for X axis
%       ytile       number of slice tiles for Y axis
%       nslices     number of slices in output volume
%       xres        X pixel resolution in mm
%       yres        Y pixel resolution in mm
%       zres        Z resolution (spacing between slice centers)
%       flip        char argument, can "contain" either of x, y, or z
%       dtype       input data type (either char or empty array of type)
%       odtype      char name of datatype
%
% Note: all arguments are mandatory!
%
% See also spm_create_vol, spm_write_vol.

% Version:  v0.5c
% Build:    6120415
% Date:     Dec-04 2006, 3:15 PM CET
% Author:   Jochen Weber, Brain Innovation, B.V., Maastricht, NL
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 13 || ...
   ~ischar(rawfile) || ...
    isempty(rawfile) || ...
    exist(rawfile, 'file') ~= 2 || ...
   ~ischar(volfile) || ...
    isempty(volfile) || ...
    strcmp(rawfile, volfile)
    error( ...
        'BVQXtools:BadArgument', ...
        'Too few arguments, invalid rawfile or volfile argument.' ...
    );
end
if ~isa(xpix, 'double') || ...
    numel(xpix) ~= 1 || ...
    xpix ~= fix(xpix) || ...
   ~isa(ypix, 'double') || ...
    numel(ypix) ~= 1 || ...
    ypix ~= fix(ypix)
    error( ...
        'BVQXtools:BadArgument', ...
        'Invalid xpix or ypix argument.' ...
    );
end
numpix = xpix * ypix;
if numpix < 1 || ...
    isnan(numpix) || isinf(numpix)
    error( ...
        'BVQXtools:BadArgument', ...
        'Invalid xpix or ypix argument.' ...
    );
end
if ~isa(xtile, 'double') || ...
    numel(xtile) ~= 1 || ...
    xtile ~= fix(xtile) || ...
   ~isa(ytile, 'double') || ...
    numel(ytile) ~= 1 || ...
    ytile ~= fix(ytile)
    error( ...
        'BVQXtools:BadArgument', ...
        'Invalid xtile or ytile argument.' ...
    );
end
numtiles = xtile * ytile;
if numtiles < 1 || ...
    isnan(numtiles) || ...
    isinf(numtiles)
    error( ...
        'BVQXtools:BadArgument', ...
        'Invalid xtile or ytile argument.' ...
    );
end
xpixs = xpix / xtile;
ypixs = ypix / ytile;
if xpixs ~= fix(xpixs) || ...
    ypixs ~= fix(ypixs)
    error( ...
        'BVQXtools:BadArgument', ...
        'Bad tiling: xtile/ytile must be a divisor of xpix/ypix.' ...
    );
end
if ~isa(nslices, 'double') || ...
    numel(nslices) ~= 1 || ...
    nslices ~= fix(nslices) || ...
    isnan(nslices) || ...
    isinf(nslices) || ...
    nslices > numtiles
    error( ...
        'BVQXtools:BadArgument', ...
        'Invalid number of slices.' ...
    );
end
if ~isa(xres, 'double') || ...
    numel(xres) ~= 1 || ...
    isnan(xres) || ...
    isinf(xres) || ...
   ~isa(yres, 'double') || ...
    numel(yres) ~= 1 || ...
    isnan(yres) || ...
    isinf(yres) || ...
   ~isa(zres, 'double') || ...
    numel(zres) ~= 1 || ...
    isnan(zres) || ...
    isinf(zres)
    error( ...
        'BVQXtools:BadArgument', ...
        'Invalid xres, yres, or zres argument.' ...
    );
end
flipx = 1;
flipy = 1;
flipz = 1;
if ischar(flip)
    flip = lower(flip(:)');
    if any(flip == 'x')
        flipx = -1;
    end
    if any(flip == 'y')
        flipy = -1;
    end
    if any(flip == 'z')
        flipz = -1;
    end
end
if ~ischar(dtype) || ...
    isempty(dtype)
    dtype = class(dtype);
end
if ~ischar(odtype) || ...
    isempty(odtype)
    odtype = class(odtype);
end

% get global defaults
global defaults;
if isempty(defaults) || ...
   ~isstruct(defaults) || ...
   ~isfield(defaults, 'analyze')
    spm_defaults;
end
try
    if defaults.analyze.flip
        flipx = flipx * -1;
    end
catch
    % nothing
end

% load pixel data
try
    vf = fopen(rawfile(:)', 'rb');
    frewind(vf);
    vd = double(fread(vf, [xpix, ypix], ['*' dtype]));
catch
    error( ...
        'BVQXtools:ErrorReadingFile', ...
        'Error reading raw data image file.' ...
    );
end

% check dimensions
if ndims(vd) > 2 || ...
    any(size(vd) ~= [xpix, ypix])
    error( ...
        'BVQXtools:DimensionMismatch', ...
        'The data read from the file doesn''t match xpix/ypix dimensions.' ...
    );
end

% build output image
try
    V = struct;
    V.fname = volfile(:)';
    V.dim   = [xpixs, ypixs, nslices, spm_type(odtype)];
    V.mat   = eye(4);
    V.mat(1,1) = xres * flipx;
    V.mat(2,2) = yres * flipy;
    V.mat(3,3) = zres * flipz;
    V.mat(1,4) = (-(xpixs   + 1) / 2) * V.mat(1,1);
    V.mat(2,4) = (-(ypixs   + 1) / 2) * V.mat(2,2);
    V.mat(3,4) = (-(nslices + 1) / 2) * V.mat(3,3);
    V.pinfo = [1;0;0];
    V = spm_create_vol(V);
    Vd = zeros([xpixs, ypixs, nslices]);
catch
    error( ...
        'BVQXtools:SPMError', ...
        'Error calling spm_create_vol for ''%s''.', ...
        V.fname ...
    );
end

% fill output image
sc = 1;
for xtc = 1:xpixs:xpix
    for ytc = 1:ypixs:ypix
        Vd(:,:,sc) = vd(ytc:(ytc+ypixs-1),xtc:(xtc+xpixs-1));
        sc = sc + 1;
        if sc > nslices, break; end
    end
    if sc > nslices, break; end
end

% write output image
try
    V = spm_write_vol(V, Vd);
    V = spm_close_vol(V);
    fclose('all');
catch
    error( ...
        'BVQXtools:SPMError', ...
        'Error calling spm_write_vol / spm_close_vol for ''%s''.', ...
        V.fname ...
    );
end
